# Clear all
rm(list = ls()) 
set.seed(47474747)

#Load Packages
library(lpSolveAPI)
library(readstata13)
library(dplyr)
library(ggplot2)
library(AER)
library(sandwich)
library(lmtest)
library(boot)

# Set Path <-- Direct to replication data folder
setwd("/Users/vedantvohra/Library/CloudStorage/Dropbox/ACR Bounds/Replication Kit")

# Load Function "cce"
source("Functions/cce.R")

# The function takes arguments (concavity_assumption, weights, acr).
# 1. concavity_assumption = 1 if beta_{m} > beta_{m+1}, 0 otherwise
# 2. weights is the vector type containing the appropriate weights
# 3. acr should be the appropriate value type containing the ACR
# Outputs Lower Bound, Upper Bound, and values of betas in both cases

# Input observed weights as a vector
ohie_weights_all <- read.dta13("Output/ohie_ed_w_all.dta") 
ohie_weights_all <- as.vector(ohie_weights_all[,2])

ohie_weights_zero <- read.dta13("Output/ohie_ed_w_zero.dta") 
ohie_weights_zero <- as.vector(ohie_weights_zero[,2])

ohie_weights_nonzero <- read.dta13("Output/ohie_ed_w_nonzero.dta") 
ohie_weights_nonzero <- as.vector(ohie_weights_nonzero[,2])

# Specify ACRs (absolute value)
acr_all <- 0.52931 
acr_zero <- 0.51662 
acr_nonzero <- 0.5565

n_all <- 24646
n_zero <- 16930
n_nonzero <- 7716 

# LP Result assuming concavity
results_conc_zero <- cce(1,ohie_weights_zero,acr_zero)
results_conc_nonzero <- cce(1,ohie_weights_nonzero,acr_nonzero)

cce_lb_all <- (n_nonzero*results_conc_nonzero$`Lower Bound` + n_zero*results_conc_zero$`Lower Bound` )/n_all
cce_ub_all <- (n_nonzero*results_conc_nonzero$`Upper Bound` + n_zero*results_conc_zero$`Upper Bound` )/n_all

source("Functions/cce_j1_j2.R")

# LP Result for 12 months
results_conc_zero <- cce_j1_j2(1,ohie_weights_zero,acr_zero,0,12)
results_conc_nonzero <- cce_j1_j2(1,ohie_weights_nonzero,acr_nonzero,0,12)

#12 Months
cce_lb_all_12 <- (n_nonzero*results_conc_nonzero$`Lower Bound` + n_zero*results_conc_zero$`Lower Bound` )/n_all
cce_ub_all_12 <- (n_nonzero*results_conc_nonzero$`Upper Bound` + n_zero*results_conc_zero$`Upper Bound` )/n_all

ub<-c()
lb<-c()
dosage<-c()

for (i in 1:19){
  
  # LP Result for i months
  results_conc_zero <- cce_j1_j2(1,ohie_weights_zero,acr_zero,0,i)
  results_conc_nonzero <- cce_j1_j2(1,ohie_weights_nonzero,acr_nonzero,0,i)
  
  # Combine for zero and non-zero ED visits in the past
  cce_lb <- (n_nonzero*results_conc_nonzero$`Lower Bound` + n_zero*results_conc_zero$`Lower Bound` )/n_all
  cce_ub <- (n_nonzero*results_conc_nonzero$`Upper Bound` + n_zero*results_conc_zero$`Upper Bound` )/n_all
  
  lb<- c(lb,cce_lb)
  ub<- c(ub,cce_ub)
  
  dosage <- c(dosage,i)
}

data <- data_frame(ub,lb,dosage)

plot <- ggplot(data,aes(x = dosage)) +
  geom_line(aes(y = ub, color = "Upper Bound")) +
  geom_line(aes(y = lb, color = "Lower Bound")) +
  labs(y=expression("CE (0, "*j[2]*")"),x =expression("Months of Medicaid ("*j[2]*")")) +
  scale_color_manual(values = c("blue", "red"), name = NULL) + 
  theme_minimal() +
  scale_x_continuous(breaks = seq(1, 19, by = 1)) +
  theme(panel.grid.minor = element_blank()) 
plot
ggsave("Output/ohie_cce_varying_j.jpg", plot, width = 8, height = 6, dpi = 300)

